function tables_compstat_SES(col_long,fp,output_dir)

% Comparative statics for SES-dependent parameters:
% 1) baseline 1: full heterogeneity
% 2) baseline 2: shut down age-dependence in wages (but everything else varies by SES)
% 3) assume homogeneous wages
% 4) homogeneous wages and tech
% 5) homogeneous wages, tech and discount factors

% Load appropriate results
load([output_dir 'temp_compstat_SES'],'sim_data_all');

% Number of models compared:
nr_models = size(sim_data_all,3);
% Slice 4 = baseline

% Tables containing all the averages, for the baseline cases and the comparative statics cases
table_means = NaN(28,nr_models); % means across all ages
table_gaps = NaN(28,nr_models); % absolute and percentage gaps
tmp_base = 1; % define baseline model for percentage gaps
% 1 = baseline model with full heterogeneity,
% 2 = baseline model where parental wages no longer vary with parental age

% Define cutoffs for comparing "low SES" to "high SES"
educ_low = 12;  % less than or equal to this
educ_high = 17; % more than or equal to this

for i = 1:nr_models
    % define baselines for each of the models (without policy intervention, so the ''end'')
    tmp_data = sim_data_all(:,:,i);
    
    % Define proper rows for period M+1 outcomes
    tmprowsL = tmp_data(:,col_long.age) == fp.M & tmp_data(:,col_long.dad_educ) <= educ_low;
    tmprowsH = tmp_data(:,col_long.age) == fp.M & tmp_data(:,col_long.dad_educ) >= educ_high;
    
    % 1) and 2) log k_17, mean for low vs high SES
    table_means(1,i) = nanmean(log(tmp_data(tmprowsL, col_long.k_tp1))); % latent child capital at age fp.M+1 = 17 (off grid values)
    table_means(2,i) = nanmean(log(tmp_data(tmprowsH, col_long.k_tp1))); % latent child capital at age fp.M+1 = 17 (off grid values)
    
    table_gaps(1,i) = table_means(2,i)-table_means(1,i); % absolute gap: high SES - low SES
    table_gaps(2,i) = 100*table_gaps(1,i) / table_gaps(1,tmp_base); % percentage gap: (high SES - low SES) / baseline gap
    
    % 3) and 4) final beta_c,17, low vs high SES
    tmp0 = tmp_data(tmprowsL, col_long.n_tp1); % latent child capital at age fp.M+1 = 17 (off grid values)
    table_means(3,i) = nanmean(tmp0./(1+tmp0));% transform n into beta_c
    
    tmp0 = tmp_data(tmprowsH, col_long.n_tp1); % latent child capital at age fp.M+1 = 17 (off grid values)
    table_means(4,i) = nanmean(tmp0./(1+tmp0)); % transform n into beta_c
    
    table_gaps(3,i) = table_means(4,i)-table_means(3,i); % absolute gap: high SES - low SES
    table_gaps(4,i) = 100*table_gaps(3,i) / table_gaps(3,tmp_base); % percentage gap: (high SES - low SES) / baseline gap
    
    % Define new rows for ages 3-16 of other outcomes/choices
    tmprowsL = tmp_data(:,col_long.age) >= 3 & tmp_data(:,col_long.age) <= 16 & tmp_data(:,col_long.dad_educ) <= educ_low;
    tmprowsH = tmp_data(:,col_long.age) >= 3 & tmp_data(:,col_long.age) <= 16 & tmp_data(:,col_long.dad_educ) >= educ_high;
    
    % 5) and 6) CCT use, low vs high SES
    table_means(5,i) = nanmean(tmp_data(tmprowsL, col_long.CCT));
    table_means(6,i) = nanmean(tmp_data(tmprowsH, col_long.CCT));
    
    table_gaps(5,i) = table_means(6,i)-table_means(5,i); % absolute gap: high SES - low SES
    table_gaps(6,i) = 100*table_gaps(5,i) / table_gaps(5,tmp_base); % percentage gap: (high SES - low SES) / baseline gap
    
    % 13) and 14) w1, low vs high SES
    table_means(13,i) = nanmean(tmp_data(tmprowsL, col_long.mom_wage));
    table_means(14,i) = nanmean(tmp_data(tmprowsH, col_long.mom_wage));
    
    table_gaps(13,i) = table_means(14,i)-table_means(13,i); % absolute gap: high SES - low SES
    table_gaps(14,i) = 100*table_gaps(13,i) / table_gaps(13,tmp_base); % percentage gap: (high SES - low SES) / baseline gap
    
    % 15) and 16) w2, low vs high SES
    table_means(15,i) = nanmean(tmp_data(tmprowsL, col_long.dad_wage));
    table_means(16,i) = nanmean(tmp_data(tmprowsH, col_long.dad_wage));
    
    table_gaps(15,i) = table_means(16,i)-table_means(15,i); % absolute gap: high SES - low SES
    table_gaps(16,i) = 100*table_gaps(15,i) / table_gaps(15,tmp_base); % percentage gap: (high SES - low SES) / baseline gap

    % 27) and 28) tau_p = tau1+tau2, low vs high SES
    table_means(27,i) = nanmean(tmp_data(tmprowsL, col_long.tau1)+tmp_data(tmprowsL, col_long.tau2));
    table_means(28,i) = nanmean(tmp_data(tmprowsH, col_long.tau1)+tmp_data(tmprowsH, col_long.tau2));
    
    table_gaps(27,i) = table_means(28,i)-table_means(27,i); % absolute gap: high SES - low SES
    table_gaps(28,i) = 100*table_gaps(27,i) / table_gaps(27,tmp_base); % percentage gap: (high SES - low SES) / baseline gap
end

%% Now create a table with the summary statistics for each alternative model + baseline in rightmost column

tmp = [table_means(:,1) table_means(:,2) table_means(:,3) table_means(:,4) table_means(:,5)];
tmp2 = [table_gaps(:,1) table_gaps(:,2) table_gaps(:,3) table_gaps(:,4) table_gaps(:,5)];

diary off;
diary_dir = [cd,'\table_compstat_SES.txt'];
force_delete(diary_dir);
diary('table_compstat_SES.txt');
diary on;

prec = '%.3f'; % precision for this table
prec2 = '%.1f'; % precision for this table

disp(['\begin{table}[!htb]']);
disp(['\begin{center}']);
disp(['\caption{Comparative Statics - Heterogeneity based on Demographics}']);
disp(['\label{table:compstat_SES}']);
disp(['\begin{tabular*}{\textwidth}{@{\extracolsep{\fill}}lccccc}']);
disp(['\hline \hline']);
disp([' & \multicolumn{1}{c}{\textbf{Baseline}} \rule{0pt}{3ex}',...
    '& \multicolumn{1}{c}{\textbf{(1)}} ',...
    '& \multicolumn{1}{c}{\textbf{(2)}} ',...
    '& \multicolumn{1}{c}{\textbf{(3)}} ',...
    '& \multicolumn{1}{c}{\textbf{(4)}} \\']);
disp(['\hline']);
disp([' \textbf{Model specifications} \rule{0pt}{3ex} & & & & \\']);
disp(['\quad - Wages/NLI = f(age)\rule{0pt}{3ex}          & Yes & No  & No    & No    & No \\']);
disp(['\quad - Wages/NLI = f(educ.)        & Yes & Yes & No    & No    & No \\']);
disp(['\quad - Time prod. = f(educ.)        & Yes & Yes & Yes   & No    & No \\']);
disp(['\quad - Discount factor distr. = f(educ.)   & Yes & Yes & Yes   & Yes   & No \\']);
disp(['\quad - Fix child''s initial test score   & Yes & Yes & Yes   & Yes   & Yes \\']);

disp(['\hline']);
disp([' \textbf{Final-Period Skills} \rule{0pt}{3ex} & & & & \\']);
disp(['Final-period Child Quality $\log(k_{17})$, Low SES \rule{0pt}{3ex}& ',num2str(tmp(1,1),prec),' & ',num2str(tmp(1,2),prec),' & ',num2str(tmp(1,3),prec),' & ',num2str(tmp(1,4),prec),' & ',num2str(tmp(1,5),prec),' \\ ', ...
    'Final-period Child Quality $\log(k_{17})$, High SES & ',num2str(tmp(2,1),prec),' & ',num2str(tmp(2,2),prec),' & ',num2str(tmp(2,3),prec),' & ',num2str(tmp(2,4),prec),' & ',num2str(tmp(2,5),prec),' \\ ', ...
    'Gap relative to baseline (\%) & ',num2str(tmp2(2,1),prec2),' & ',num2str(tmp2(2,2),prec2),' & ',num2str(tmp2(2,3),prec2),' & ',num2str(tmp2(2,4),prec2),' & ',num2str(tmp2(2,5),prec2),' \\ ', ...
    'Final-period Discount Factor $\beta_{c,17}$, Low SES \rule{0pt}{3ex}& ',num2str(tmp(3,1),prec),' & ',num2str(tmp(3,2),prec),' & ',num2str(tmp(3,3),prec),' & ',num2str(tmp(3,4),prec),' & ',num2str(tmp(3,5),prec),' \\ ',...
    'Final-period Discount Factor $\beta_{c,17}$, High SES & ',num2str(tmp(4,1),prec),' & ',num2str(tmp(4,2),prec),' & ',num2str(tmp(4,3),prec),' & ',num2str(tmp(4,4),prec),' & ',num2str(tmp(4,5),prec),' \\ ',...
    'Gap relative to baseline (\%)  & ',num2str(tmp2(4,1),prec2),' & ',num2str(tmp2(4,2),prec2),' & ',num2str(tmp2(4,3),prec2),' & ',num2str(tmp2(4,4),prec2),' & ',num2str(tmp2(4,5),prec2),' \\ ']);

disp([' \textbf{Parental Choices} \rule{0pt}{3ex} & & & & \\']);

disp(['Fraction using CCT, Low SES \rule{0pt}{3ex} & ',num2str(tmp(5,1),prec),' & ',num2str(tmp(5,2),prec),' & ',num2str(tmp(5,3),prec),' & ',num2str(tmp(5,4),prec),' & ',num2str(tmp(5,5),prec),' \\ ',...
    'Fraction using CCT, High SES & ',num2str(tmp(6,1),prec),' & ',num2str(tmp(6,2),prec),' & ',num2str(tmp(6,3),prec),' & ',num2str(tmp(6,4),prec),' & ',num2str(tmp(6,5),prec),' \\ ',...
    'Gap relative to baseline (\%)  & ',num2str(tmp2(6,1),prec2),' & ',num2str(tmp2(6,2),prec2),' & ',num2str(tmp2(6,3),prec2),' & ',num2str(tmp2(6,4),prec2),' & ',num2str(tmp2(6,5),prec2),' \\ ',...
    'Parental investment time ($\tau_p$), Low SES \rule{0pt}{3ex} & ',num2str(tmp(27,1),prec),' & ',num2str(tmp(27,2),prec),' & ',num2str(tmp(27,3),prec),' & ',num2str(tmp(27,4),prec),' & ',num2str(tmp(27,5),prec),' \\ ',...
    'Parental investment time ($\tau_p$), High SES & ',num2str(tmp(28,1),prec),' & ',num2str(tmp(28,2),prec),' & ',num2str(tmp(28,3),prec),' & ',num2str(tmp(28,4),prec),' & ',num2str(tmp(28,5),prec),' \\ ',...
    'Gap relative to baseline (\%)  & ',num2str(tmp2(28,1),prec2),' & ',num2str(tmp2(28,2),prec2),' & ',num2str(tmp2(28,3),prec2),' & ',num2str(tmp2(28,4),prec2),' & ',num2str(tmp2(28,5),prec2),' \\ ']);

disp([' \textbf{Parental Wages} \rule{0pt}{3ex} & & & & \\']);

disp(['Mother''s hourly wage, Low SES \rule{0pt}{3ex} & ',num2str(tmp(13,1),prec),' & ',num2str(tmp(13,2),prec),' & ',num2str(tmp(13,3),prec),' & ',num2str(tmp(13,4),prec),' & ',num2str(tmp(13,5),prec),' \\ ',...
    'Mother''s hourly wage, High SES & ',num2str(tmp(14,1),prec),' & ',num2str(tmp(14,2),prec),' & ',num2str(tmp(14,3),prec),' & ',num2str(tmp(14,4),prec),' & ',num2str(tmp(14,5),prec),' \\ ',...
    'Gap relative to baseline (\%)  & ',num2str(tmp2(14,1),prec2),' & ',num2str(tmp2(14,2),prec2),' & ',num2str(tmp2(14,3),prec2),' & ',num2str(tmp2(14,4),prec2),' & ',num2str(tmp2(14,5),prec2),' \\ ',...
    'Father''s hourly wage, Low SES \rule{0pt}{3ex} & ',num2str(tmp(15,1),prec),' & ',num2str(tmp(15,2),prec),' & ',num2str(tmp(15,3),prec),' & ',num2str(tmp(15,4),prec),' & ',num2str(tmp(15,5),prec),' \\ ',...
    'Father''s hourly wage, High SES & ',num2str(tmp(16,1),prec),' & ',num2str(tmp(16,2),prec),' & ',num2str(tmp(16,3),prec),' & ',num2str(tmp(16,4),prec),' & ',num2str(tmp(16,5),prec),' \\ ',...
    'Gap relative to baseline (\%)  & ',num2str(tmp2(16,1),prec2),' & ',num2str(tmp2(16,2),prec2),' & ',num2str(tmp2(16,3),prec2),' & ',num2str(tmp2(16,4),prec2),' & ',num2str(tmp2(16,5),prec2),' \\ ',...
    ]);

disp(['\hline \hline']);
disp(['\end{tabular*}']);
disp(['\end{center}']);
disp(['\break']);
disp(['\textbf{Notes:} ``Baseline" corresponds to the baseline model. ', ... % where all primitives (i.e., parental wages, household discount factors and parental time productivities) vary with parental education. ',...Column (1) corresponds to the alternative model where all model primitives are assumed to be independent of parental education. \\',...
    'Column (1) corresponds to the alternative model where parental wage offers and non-labor income (NLI) no longer depend on parental age. ',...
    'Column (2) considers the case where parental wages and NLI no longer depend on parental age or education levels. ',...
    'Column (3) considers the case where parental wages, NLI and investment time productivities no longer depend on parental characteristics. ',...
    'Column (4) considers the case where parental wages, NLI, investment time productivities, and parent/child (initial) discount factor distributions no longer depend on parental characteristics. ',...
    '``Low SES" considers simulated households where the father has at most a high school degree, whereas ``High SES" looks at households where the father has a graduate degree. ',...
    'All experiments were done using $R = ',num2str(fp.R),'$ simulated data sets.']);
disp(['\end{table}']);

diary off;

end % function